{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "function [L,U] = LU(A)\n",
    "    n = length(A);\n",
    "    L = eye(n);\n",
    "    U = zeros(n);\n",
    "    for r = 1:n\n",
    "        U(r,r:n) = A(r,r:n) - L(r,1:r-1)*U(1:r-1,r:n);\n",
    "        L(r+1:n,r) = (A(r+1:n,r) - L(r+1:n,1:r-1)*U(1:r-1,r))/U(r,r);\n",
    "    end\n",
    "endfunction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =\n",
      "\n",
      "   10.00000   -7.00000    0.00000    1.00000\n",
      "   -3.00000    2.10000    6.00000    2.00000\n",
      "    5.00000   -1.00000    5.00000   -1.00000\n",
      "    2.00000    1.00000    0.00000    2.00000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "A = [10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans = [](1x0)\n",
      "ans = [](0x4)\n",
      "ans = -0.30000\n",
      "ans =\n",
      "\n",
      "  -7   0   1\n",
      "\n",
      "ans =\n",
      "\n",
      "         0.50000  -2499999.99965\n",
      "\n",
      "ans =\n",
      "\n",
      "   0.00000   1.00000\n",
      "   6.00000   2.30000\n",
      "\n",
      "ans =\n",
      "\n",
      "         0.20000  -2399999.99966         0.96000\n",
      "\n",
      "ans =\n",
      "\n",
      "         1.00000\n",
      "         2.30000\n",
      "   5749998.49920\n",
      "\n"
     ]
    }
   ],
   "source": [
    "[L,U]=LU(A);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L =\n",
      "\n",
      "         1.00000         0.00000         0.00000         0.00000\n",
      "        -0.30000         1.00000         0.00000         0.00000\n",
      "         0.50000  -2499999.99965         1.00000         0.00000\n",
      "         0.20000  -2399999.99966         0.96000         1.00000\n",
      "\n",
      "U =\n",
      "\n",
      "         10.00000         -7.00000          0.00000          1.00000\n",
      "          0.00000         -0.00000          6.00000          2.30000\n",
      "          0.00000          0.00000   15000004.99790    5749998.49920\n",
      "          0.00000          0.00000          0.00000          5.08000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "L,U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "function x = LinearSolveWithLU(L,U,b)\n",
    "    n = length(b);\n",
    "    y = zeros(n,1);\n",
    "    for i = 1:n   \n",
    "        y(i) = b(i)-L(i,1:i-1)*y(1:i-1);\n",
    "    end\n",
    "    x = zeros(n,1);\n",
    "    for i = n:-1:1\n",
    "        x(i) = (y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);\n",
    "    end\n",
    "endfunction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b =\n",
      "\n",
      "   8.0000\n",
      "   5.9000\n",
      "   5.0000\n",
      "   1.0000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "b = [8;5.900001;5;1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x =\n",
      "\n",
      "  -0.00000000060339\n",
      "  -1.00000000088818\n",
      "   1.00000000007028\n",
      "   0.99999999981667\n",
      "\n"
     ]
    }
   ],
   "source": [
    "x = LinearSolveWithLU(L,U,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "5.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
